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Abstract 

New  methods  for  the  analysis  of  three-component  seismograms  have  been  applied  to  data  from 
both  continental  and  oceanic  regions.  Polarization  anisotropy,  manifested  as  the  splitting  of 
surface,  guided,  and  shear  body  waves,  has  been  observed  in  all  regions,  but  no  significant 
azimuthal  anisotropy  has  been  detected.  This  is  expected  if  the  local  orientation  of  olivine 
crystals  in  the  upper  mantle  is  incoherent  on  the  scale  of  the  path  lengths  used  in  these 
experiments  (A  =  35°-50°).  From  a  detailed  study  of  trans-Australia  paths,  we  have  shown  that 
the  shear-wave  anisotropy  in  the  western  Australia  craton  averages  3-4%  in  the  uppermost 
mantle,  but  it  terminates  abruptly  at  a  Lehmann  (L)  discontinuity  near  a  depth  of  250  km.  Our 
physical  interpretation  is  that  the  L  discontinuity  beneath  the  ancient  continental  cratons  marks 
the  transition  from  an  anisotropic  mechanical  boundary  layer  to  a  more  isotropic  region  of  a  thick 
continental  tectosphere.  Above  L,  the  temperatures  have  evidently  remained  cold  enough  to 
freeze-in  the  small-scale  anisotropic  structures  and  tectonic  fabrics  generated  in  episodes  of 
orogenic  compression  associated  with  tectospheric  stabilization.  Below  L,  such  structures  were 
either  never  generated  or  were  annealed  out  subsequent  to  their  formation.  These  hypotheses 
about  subcontinental  structure  have  considerable  bearing  on  a  diverse  set  of  geological  and 
geodynamical  problems,  ranging  from  formation  of  the  cratons  to  kimberlite  vulcanism. 

To  provide  additional  quantitative  tests  of  these  hypotheses,  we  have  formulated  stochastic 
models  that  specify  the  small-scale  heterogeneity  of  the  continental  upper  mantle  as  samples  of 
Gaussian  random  fields  with  self-affine  (fractal)  scaling  at  high  wavenumbers.  We  approximate 
the  spatial  variations  in  elasticity  C(x)  as  a  hexagonal  tensor  field  with  a  local  axis  of  symmetry 
given  by  a  unit  vector  field  s(x).  We  impose  the  structure  of  a  stationary  Gaussian  random  field 
on  s(x)  and,  from  a  deterministic  functional  relationship  between  C  and  s,  calculate  the 
ensemble  averages  of  C(x)  that  are  needed  to  describe  low-frequency  wave  propagation. 
Solving  forward  problems  of  this  type  allows  us  to  set  up  quantitative  inverse  problems  for  the 
parameters  of  the  stochastic  distributions.  In  the  modeling  discussed  here,  s  is  assumed  to  have 
transversely  isotropic  stochastic  symmetry  and  to  be  specified  by  four  statistical  parameters:  an 
aspect  ratio  of  the  anisotropy,  f,  a  characteristic  horizontal  wavenumber  of  the  heterogeneity,  k, 
an  aspect  ratio  of  the  heterogeneity,  rj,  and  a  fractal  dimension,  D.  For  =  1 ,  the  distribution  of 
s  is  isotropic  in  all  three  directions;  in  the  limit  §  — »  °o ,  s3  equals  zero  with  probability  one,  and 
s  is  isotropically  distributed  in  the  horizontal  plane.  We  show  that  the  upper-mantle  anisotropy 
in  Gaherty  and  Jordan's  [1995]  model  AU3  is  consistent  with  Voigt  averages  given  by  f  «  2-3. 
We  discuss  how  better  approximations  to  the  effective  elastic  tensor  can  be  obtained  from  a  self- 
consistent,  second-order  Bom  approximation. 
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Objectives 


The  technologies  of  nuclear-monitoring  seismology  are  built,  in  part,  on  studies  of  wave 
propagation  at  regional  distances  (<  2000  km).  One  poorly  understood  aspect  of  regional 
propagation  is  the  forward  scattering  of  high-frequency  signals  by  small  inhomogeneities.  The 
objectives  of  this  research  are  threefold:  (1)  to  formulate  parameterized  stochastic  models  of 
small-scale  heterogeneity  that  are  suitable  for  assessing  scattering  effects  at  regional  distances, 
(2)  to  constrain  the  model  parameters  using  broadband  seismic  data,  and  (3)  to  assess  the 
implications  of  these  models  for  regional  wave  propagation  related  to  nuclear  monitoring. 

Research  Accomplishments:  Observations  and  Modeling  of  Seismic  Corridors 

Our  approach  is  to  model  small-scale  heterogeneity  by  a  bootstrap  from  low  to  high  frequencies 
[Jordan  and  Gaherty,  1994].  In  the  first  phase  of  this  program,  we  have  used  low-frequency 
(<  .05  Hz),  phase-coherent  waves  to  derive  anisotropic,  path-averaged  models  of  crustal  and 
upper-mantle  structure.  We  have  developed  methods  for  extracting  phase  delays  and  amplitudes 
from  three-component  seismograms  that  have  delivered  new  information  about  impedance 
discontinuities  and  velocity  gradients  in  the  upper  mantle  [Revenaugh  and  Jordan,  1991;  Gee  and 
Jordan,  1992],  and  we  have  applied  them  to  data  from  both  continental  [Gaherty  and  Jordan, 
1995]  and  oceanic  [Gaherty  et  al.,  1995]  regions.  Polarization  anisotropy,  manifested  as  the 
splitting  of  surface,  guided,  and  shear  body  waves,  has  been  observed  in  all  regions,  but  no 
significant  azimuthal  anisotropy  has  been  detected.  This  is  expected  if  the  local  orientation  of 
olivine  crystals  in  the  upper  mantle  is  incoherent  on  the  scale  of  the  path  lengths  used  in  these 
studies  (A  =  35°-50°).  We  have  therefore  modeled  the  splitting  in  terms  of  radially  anisotropic, 
but  transversely  isotropic,  elasticity.  From  a  detailed  study  of  trans-Australia  paths  [Gaherty  and 
Jordan,  1995],  we  have  shown  that  the  shear-wave  anisotropy  in  the  western  Australia  craton 
averages  3-4%  in  the  uppermost  mantle,  but  it  terminates  abruptly  at  a  Lehmann  (L) 
discontinuity  near  a  depth  of  250  km.  At  the  L  transition,  which  occurs  over  a  depth  range  of 
less  than  30  km,  the  average  SV  wavespeed  increases  and  the  average  SH  wavespeed  decreases, 
obtaining  nearly  equal  values  in  the  mantle  layer  between  L  and  the  410-km  discontinuity. 

Our  physical  interpretation  is  that  the  L  discontinuity  beneath  the  ancient  continental  cratons 
marks  the  transition  from  an  anisotropic  mechanical  boundary  layer  to  a  more  isotropic  region  of 
a  thick  continental  tectosphere  [Revenaugh  and  Jordan,  1991;  Gaherty  and  Jordan,  1995].  Above 
L,  the  temperatures  have  evidently  remained  cold  enough  to  freeze-in  the  small-scale  anisotropic 
structures  and  tectonic  fabrics  generated  in  episodes  of  orogenic  compression  associated  with 
tectospheric  stabilization  [Jordan,  1988;  Silver  and  Chan,  1991],  Below  L,  such  structures  were 
either  never  generated,  perhaps  because  L  coincides  with  a  transition  from  dislocation  to 
diffusion  creep  [Karato,  1992],  or  were  annealed  out  subsequent  to  their  formation.  These 
hypotheses  about  subcontinental  structure  have  considerable  bearing  on  a  diverse  set  of 
geological  and  geodynamical  problems,  ranging  from  formation  of  the  cratons  to  kimberlite 
vulcanism  [Gaherty  and  Jordan,  1995]. 


434 


Research  Accomplishments:  Stochastic  Modeling  of  Fine-Scale  Anisotropic  Structures 

In  the  current  phase  of  our  research  program,  we  are  focusing  on  the  interpretation  of 
polarization  anisotropy  in  terms  of  heterogeneity  at  subwavelength  scales.  We  have  formulated 
stochastic  models  that  specify  the  small-scale  heterogeneity  of  the  continental  upper  mantle  as 
samples  of  Gaussian  random  fields  with  self-affine  (fractal)  scaling  at  high  wavenumbers,  and 
we  have  derived  expressions  relating  the  parameters  of  these  models  to  the  speeds  of  low- 
frequency  seismic  waves.  This  report  presents  some  of  these  results  and  outlines  future  research 
activity  related  to  this  modeling. 

We  consider  a  medium  governed  by  a  linear  stress-strain  relation  otJ  =  Cl]kiEkl ,  where  the 
elasticity  tensor  C  has  Cartesian  components  {C^fc/x) :  i,j,k,l  =  1,2,3}  that  are  functions  of 
three-dimensional  position  vector  x.  Theoretical  calculations  [Estey  and  Douglas,  1986]  and 
measurements  on  kimberlite  xenoliths  [Mainprice  and  Silver,  1993;  Christensen,  1994]  indicate 
that  the  elasticity  of  upper-mantle  peridotites  has  quasi-hexagonal  symmetry  on  the  hand-sample 
scale.  Therefore,  it  is  reasonable  to  approximate  C(x)  as  a  hexagonal  tensor  field  with  a  local 
axis  of  symmetry  given  by  the  unit  vector  s(x)  =  s(x)/ls(x)l.  This  can  be  accomplished  by 
restricting  the  tensor  field  to  be  a  functional  of  a  vector  field  s(x):  C(x)  =  C(s(x)).  In 
petrological  terms,  s(x)  is  largely  determined  by  the  local  average  orientation  of  the  a  axis  of 
olivine,  which  is  aligned  by  mantle  deformation  [Estey  and  Douglas,  1986].  Our  strategy  is  to 
impose  the  structure  of  a  Gaussian  random  field  on  s(x)  and,  from  the  functional  relationship 
C(s),  to  calculate  the  ensemble  averages  of  C(x)  that  are  needed  to  describe  low-frequency 
wave  propagation.  Solving  this  forward  problem  will  allow  us  to  assess  the  applicability  of  the 
stochastic  models  to  the  continental  upper  mantle  and,  if  justified,  to  set  up  quantitative  inverse 
problems  for  the  parameters  of  the  stochastic  distributions. 

Specification  of  the  Vector  Field.  Within  some  region  of  the  upper  mantle,  the  Gaussian  random 
field  s(x)  is  assumed  to  be  stationary  with  zero  mean,  (s(x))  =0,  and  thus  specified  by  its 
autocovariance  matrix,  Css(x)  =  (s(x1  )sT(x'  +x)).  Css  is  symmetric  and  positive  definite  and 
admits  the  zero-lag  eigenvector  expansion 

3 

Cs,(0).S02-2Xe^.  (1) 

n= 1 

We  specialize  this  matrix  in  the  following  way:  its  last  eigenvector  is  assumed  to  be  vertical, 
e3  =  z,  with  dimensionless  variance  f 2 ,  and  the  variances  of  the  other  two  eigenvectors  are  set 
equal  to  unity.  Hence,  the  s  field  is  required  to  be  statistically  isotropic  in  the  xy  plane,  with 
principal  components  having  root-mean-square  (r.m.s.)  variations  given  by 

°i  -  °2  =  §  =  1  •  (2) 

We  will  call  the  dimensionless  parameter  §  the  aspect  ratio  of  the  anisotropy.  In  terms  of  the 
coordinates  5j  =  5  sin  Ocoscp ,  s2  =  ssin0sin<p,  s3  =  scosd ,  the  three-dimensional  probability 
density  element  (p.d.e.)  for  s  can  be  written 
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f  -r(*i2+s>2+i2*32) 

p(s)dV(s)  =  '^"372  e  ds1ds2ds3 


_  -s2I2  k(s)cos26  2,  ■ 

~  a  '  '  r-  /Tr  fir 


(2sc) 


3/2 


e  e“'  '  ”  '  s  dssin  ddddcp. 


K(S)  =  ~s2(  1-f2).  (3) 


From  this  expression,  it  can  be  seen  that  the  p.d.f.  on  a  sphere  of  radius  s  corresponds  to  a 
Watson  distribution  with  dispersion  parameter  k(s )  [Fisher  et  al.,  1987].  For  §  =  1,  the 
distribution  of  s  is  isotropic  in  all  three  directions.  Where  §  <  1,  the  probability  of  s  is 
concentrated  towards  the  vertical  axis  (bipolar  distribution,  k  >  0);  where  §  >  1,  it  is 
concentrated  towards  the  equatorial  plane  (girdle  distribution,  k  <  0).  In  the  limit  f  ->  o° ,  s3 
equals  zero  with  probability  one,  and  s  is  isotropically  distributed  in  the  xy  plane.  For  fixed  f, 
the  conditional  distribution  p(s  I  s)  becomes  more  isotropic  as  s  — >0. 

To  specify  the  spatial  correlation  of  the  s  field,  we  define  the  Riemannian  length. 


r2(x)  =  xtQx,  (4) 

Q  is  a  symmetric  positive-definite  matrix  chosen,  like  S0,  to  be  transversely  isotropic;  i.e.,  to 
have  eigenvectors  {e1?  e2,  e3  =  z}  and  eigenvalues 

k2  =  k\  =  rj2k2  =  k2.  (5) 

We  adopt  an  autocovariance  matrix  of  the  form 


Css(x)  =  Sq  p(r(x)) ,  (6) 

where  p(r)  is  a  scalar-valued  autocorrelation  function  such  that  p( 0)  =  1  and  lim  p(r)  =  0.  A 
convenient  and  useful  function  with  these  properties  is  r_>“ 

p(r)  =  2-(v+1)rv  £v(r)  /  T(v),  0<r<®,  0sv<l  (7) 

Here  Kv(r )  is  the  modified  Bessel  function  of  the  second  kind  and  T(  v)  is  the  gamma  function. 
This  anisotropic  autocorrelation  function  was  introduced  as  a  stochastic  representation  of 
seafloor  roughness  by  Goff  and  Jordan  [1988],  who  describe  many  of  its  properties.  A  three- 
dimensional  Fourier  transform  yields  the  power  spectrum: 

P(k)  -  f p(r(x))e~ikT x dV(x)  =  8jt3/2^^^IQI"1/2[kT Q_1k  +  1  ]"<v+3/2>  (8) 

j  r(v) 

This  spectral  density  has  principal  axes  with  comer  wavenumbers  k,  k,  t] k;  it  is  flat 

below  these  comers  and  rolls  off  above  them  at  an  asymptotic  rate  of  -(2v  +  3).  Realizations  of 
this  Gaussian  field  are  thus  characterized  by  two  outer  scales:  k~l  in  the  xy  plane,  and  ( rjk)~  in 
the  z  direction.  The  dimensionless  parameter  t]  is  termed  the  aspect  ratio  of  the  heterogeneity. 
At  scales  much  smaller  than  ( rfk)~l,  the  field  is  statistically  self-affine  with  scaling  parameter  v , 
and  it  has  a  Hausdorff  (fractal)  dimension  of 


D  =  4-  v 


(9) 
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(see  the  appendix  of  Goff  and  Jordan  [1988]  for  definitions  and  proofs).  The  special  case  v  = 
1/2  corresponds  to  the  exponential  correlation  function,  p(r)  =  eTr ,  that  has  been  commonly 
used  in  theoretical  studies  of  wave  scattering  [e.g.,  Chernov,  I960]. 

In  summary,  we  have  constructed  a  dimensionless  vector  field  s  that  has  a  transversely  isotropic 
stochastic  symmetry.  Its  statistics  are  described  by  a  single  dimensionalized  parameter,  the 
characteristic  horizontal  wavenumber  of  the  heterogeneity  k,  ([k]  =  length~),  and  three 
dimensionless  numbers:  the  aspect  ratio  of  the  anisotropy  §,  the  aspect  ratio  of  the  heterogeneity 
rj,  and  the  fractal  dimension  D.  Figures  1  and  2  illustrate  sections  taken  through  some 
realizations  of  this  Gaussian  random  field  for  k~l  =  200  km,  D  =  3.2,  and  several  values  of  §  and 
fl- 


UX)  km 

Figure  1.  A  realization  of  the  vector  field  s  for  k  1  =  200  km,  D  =  3.2,  and  fj  =  77  =  1.  The  vectors  are  calculated  on 
a  20  km  x  20  km  grid  and  displayed  as  their  projections  onto  the  sides  of  a  500  km  x  200  km  x  200  km  rectangular 
domain.  Lower  right  panel  is  a  grayshade  plot  of  the  x2  coordinate  of  s  in  the  x2-x3  plane,  showing  the  structure  of 
the  field  at  higher  resolution  than  the  vector  sampling.  In  this  case,  the  field  is  statistically  isotropic  in  all  directions 
in  both  the  vector  orientation  and  the  heterogeneity. 
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Figure  2.  Realizations  of  the  vector  field  s  for  k  -  200  km,  D  =  3.2,  and  £  =  5,  and  values  ot  the  heterogeneity 
aspect  ratio  ij  equal  to  1  (top),  5  (middle),  and  10  ( lower  ).  In  each  case,  the  left  and  nght  panels  are  plotted  with  the 
same  conventions  as  the  upper-left  and  lower-right  panels  of  Figure  1,  respectively.  The  field  in  the  top  panel  has 
statistically  isotropic  heterogeneity  but  anisotropic  v  ector  orientation.  The  field  in  the  lower  panel  corresponds  to  an 
anisotropic  stochastic  laminate  ( r]  »  1),  All  realizations  are  isotropic  in  the  x  ,-x2  plane. 


Representation  of  the  Elasticity  Field.  The  Gaussian  random  field  structure  of  s(x)  is  mapped 
onto  the  elasticity  tensor  field  C(x)  through  a  deterministic  functional  relationship  C(s).  If  C  is 
a  functional  of  s,  then  it  must  be  hexagonal  with  symmetry  axis  s.  Any  hexagonally  symmetric 
elasticity  tensor  with  components  Cljki  can  be  represented  by  six  elastic  parameters,  which  are 
taken  here  to  be 

Qlll  =  C11  =  a>  Ql22  =  c12  =  05333  =  c33  =  c’ 

Ql33  =  c13  =  /’  Q212  =  c44  =  Q313  =  c55  =  (^) 

Only  five  of  these  parameters  are  independent,  because  the  symmetry  requires 

b  =  a -2m  .  (11) 

An  alternative  and  analytically  convenient  representation  is  the  Backus  [1970]  harmonic 
decomposition: 

ClJU  -  3$ }t 4>  +  3g ',Hm  +  3 {^,Hm  +  +  4%hm.  (12) 

Here  H'  and  h  are  harmonic  polynomials  of  degree  (  ,  and  and  4jkl  are  the  f.  th  -order 
differential  operators  (  6tJ  is  the  Kronecker  delta), 

tyjkl  ~  —o',  dj  dkdh  (13a) 


%rkl  =  ~^ijdkdi  +  dkl  dt  d  s  +  6ikdj  dt  +  d  y/  <?,  dk  +  6udj  dk  + 
tyjkl  =  fyjdkl  +  dik^jl  +  dil&jk> 


$jkl  ~  ~^($ijdkdl  +  ^kt  dL  dj  ~  ~blk  dj  dt  -—djididk  ~—f>udjdk  -  —6  jkdi <?/), 


(13b) 

(13c) 

(13d) 


<4jkl  ~  &ijdkl  ~  ^ikdjl  ~  2  fyl^jk-  (^e) 

For  a  hexagonal  fourth-order  tensor  with  symmetry  axis  s ,  the  harmonic  polynomials  are 


fif(4)(r)  =  8  H4r4P4(lTr), 

(14a) 

//2)( r)  =  -2H2r2  P2(sTf), 

(14b) 

^0)(r)  =  Hq, 

(14c) 

h{2}( r)  =  -  2/22  r2  P2(sTf), 

(14d) 

-P 

11 

O 

(14e) 
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where  Pe(cosA)  is  the  Legendre  polynomial  of  order  £  and  A  =  acos(sT  f )  is  the  angle  between 
r  and  s.  As  shown  by  Backus  [1970],  the  three  polynomials  {  :  £ =0,2,4  }  represent  the 

totally  symmetric  part  of  C,  whereas  the  two  polynomials  :  £= 0,2 }  represent  its 

antisymmetric  part.  The  elastic  constants  in  (10)  are  related  to  the  coefficients  Hi  and  ht  by 


a  =  3H4  +  6H2  +  SH0 ,  (15a) 

b  =  H4  +2H2  +  H0  +  2h2  +  ho,  (15b) 

c  =  SH4-12H2+3H0,  (15c) 

/  =  -4 H4  -  H2  +  Hq  -  h,  +h0,  (15d) 

/  =  -4 H4  -H2  +  H0  -i/fe,  (15e) 

m  =  H4  +  2H2  +  Hq  -h2  -  -/b-  (15f) 


(Here  we  have  corrected  some  errors  in  equation  (34)  of  Backus  [1970].) 

To  complete  the  functional  relationship  between  C  and  s,  we  specify  the  elastic  parameters  to  be 
monomials  in  the  vector  length  s : 

He( s)  =  H£sn,  he( s)  =  hesn ,  (16) 

where  n  is  a  non-negative  integer  that  is  constant  for  all  £ .  We  call  this  additional  model 
parameter  the  strength  exponent.  In  the  case  n  =  0,  the  harmonic  coefficients  are  independent  of 
position,  and  Cijk(ix)  is  a  constant  tensor  rotated  to  have  an  axis  of  symmetry  s(x). 

Voigt  Average.  The  Voigt  average  is  the  expected  value  of  the  elasticity  tensor,  CiyW  = 
( Qjki(x )),  which  can  be  calculated  directly  from  the  expected  values  of  the  harmonic 
polynomials.  The  latter  can  be  expressed  as  sums  over  the  integrals 


Jsn  P£(cosA(s)) p(s)dV(s).  (17) 

Integrating  the  Legendre  polynomials  term-by-term  yields 

H{4) (r)  =  H4  (35 14  -  30 12  +  3 Z0" )r4  P4(zr) ,  ( 18a) 

tf(2)( r)  =  -4(3 /2"  -  Io)r2  P2(irr),  (18b) 

^(0)(r)  =  tf0/0\  (18c) 

hi2\ r)  =  -4(3/2"  -  lo)r2P2{zr),  (18d) 
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hi0\ r)  =  fiQ.  (18e) 

Hence,  the  Voigt  average  is  transversely  isotropic  with  a  vertical  symmetry  axis,  as  expected  for 
a  random  medium  whose  statistics  have  this  symmetry.  The  integrals  appearing  in  (18)  are 


/<(?)  _i!i  i+(52-i)^)'”-3>,2‘fe 

(19) 

We  give  explicit  forms  for  n  =  0,  1: 

1 0  “ 

(20a) 

(20b) 

■2cs>-»(1-3tf)- 

(20c) 

A)  -  /— — atan-Jt;  1, 

Vi  - 1 

(21a) 

h  ~  §2-ii7°~3r 

(21b) 

(21c) 


Application  to  the  Upper  Mantle.  We  have  applied  this  theory  to  the  interpretion  of  the 
anisotropy  observed  by  Gaherty  and  Jordan  [1995].  Their  model  AU3  has  a  radially  anisotropic 
region  extending  from  the  M  discontinuity  at  30  km  depth  to  the  L  discontinuity  at  252  km 
depth.  The  magnitude  of  the  anisotropy  in  such  a  transversely  isotropic  model  is  measured  by  an 
5-wave  anisotropy  ratio,  (v5W-  v sv)  /  v sv,  and  a  similar  ratio  for  P  wavespeeds.  At  the  base  of 
the  AU3  anisotropic  region,  these  ratios  are  both  about  4%.  In  comparison,  the  anisotropy  ratios 
calculated  for  a  pyrolite  mineralogy  in  which  the  olivine  and  pyroxene  crystals  have  been 
completely  aligned  by  shear  deformation  in  the  horizontal  plane  (but  randomly  oriented  within 
that  plane)  is  about  7%  [Estey  and  Douglas,  1986].  Within  the  context  of  our  stochastic  model, 
this  type  of  average  corresponds  to  §  =  °o,  describing  the  (geologically  unexpected)  situation 
where  all  the  s  vectors  are  horizontal.  As  §  decreases,  the  anisotropy  ratios  decrease,  achieving 
values  comparable  to  AU3  for  §  «  2-3.  Figure  3  shows  this  behavior  for  the  Voigt  average  in  the 
case  where  n  =  0.  An  anisotropy  aspect  ratio  of  this  magnitude  is  geologically  reasonable, 
although  it  must  be  admitted  that  not  much  is  known  from  direct  observation  about  what  f  should 
be.  The  results  in  Figure  3  at  least  suggests  that  our  approach  to  the  stochastic  modeling  of  fine- 
scale  upper-mantle  structure  might  be  on  the  right  track. 
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S-wave  anisotropy  (%) 


Figure  3.  Plot  of  the  P- wave  anisotropy,  (vpw-  vw )  /  x pv ,  versus  the  S- wave  anisotropy,  (v5H-  vJV)  /  vsv.  The 
point  with  the  1- a  error  ellipse  is  the  upper-mantle  value  from  model  AU3  of  Gaherty  and  Jordan  [1995].  The  large 
square  is  the  value  calculated  by  Voigt  averaging  the  hexagonally  symmetric  pyrolite  model  of  Estey  and  Douglas 
[1986]  assuming  complete  crystallographic  alignment  in  the  horizontal  plane  ( §  =  oo).  The  small  squares  are  Voigt 
averages  for  a  strength  exponent  n-  0  at  integer  values  of  the  anisotropy  aspect  ratio  showing  how  the  apparent 
anisotropy  decreases  as  the  local  orientation  becomes  statistically  more  isotropic.  This  calculation  suggests  that  the 
value  of  f appropriate  for  the  Australian  upper  mantle  is  on  the  order  of  2  or  3. 


Our  analysis  is  not  complete,  however.  The  one-point  probability  density  p(  s)  is  independent  of 
the  two-point  correlation  parameters;  hence,  the  Voigt  average  depends  only  on  the  aspect  ratio 
of  the  anisotropy  f  and  the  strength  exponent  n.  On  the  other  hand,  the  effective  anisotropy 
sensed  by  low-frequency  waves  must  also  be  dependent  on  the  aspect  ratio  of  the  heterogeneity 
rj.  This  illustrates  the  well  known  fact  that  the  Voigt  average  can  be  a  poor  approximation  to  the 
effective  elasticity  tensor  governing  low-frequency  wave  propagation  [Backus,  1962;  Berryman, 
1994]. 

Self-Consistent  Effective  Media  Theory:  The  Second-Order  Bom  Approximation.  Better  ap¬ 
proximations  can  be  obtained  from  self-consistent  theories  of  effective  media  based  on  an 
integral-equation  formulation  of  the  scattering  problem.  This  technique  was  developed  and  first 
applied  to  polycrystalline  aggregates  by  Korringa[1973],  and  Zeller  and  Dederichs  [1973];  see 
reviews  by  Gubematis  and  Krumhansl  [1975],  Domany  et  al.  [1975],  and  Berryman  [1994].  In 
the  long-wavelength  limit,  the  time  derivatives  of  the  dynamic  displacement  field  commute  with 
ensemble  averaging,  so  that  the  effective  elasticity  tensor 

cm  -  {o, ,)<£*,>-'  (22) 
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can  be  calculated  from  the  static  problem.  The  development  of  these  theories  begins  by 
expressing  Cljk^x)  as  the  sum  of  a  constant  tensor  C^ki  and  a  perturbation  6Cljki(x)\  then,  the 
condition  for  static  equilibrium  can  be  put  in  the  form  of  a  Lippman-Schwinger  equation 

=  ^iy(x)  +  jGijkj(x,x')dCklmn(x'  )emn(x'  )dV(x! ).  (23) 

Gl]ki  is  the  strain  Green’s  function  corresponding  to  Cijk£  x) .  In  operator  notation,  this  equation 
can  be  written  as  e  =  e°  +  G  SC  e  and  solved  iteratively  via  a  Bom  series.  The  series  can  be 
summed  to  obtain 

e  =  £°+GTe°?  (24) 

T  =  6C(I  -G&Cf1  =  6C  +  6CG6C  +  6CG6CG&C  +  •••.  (25) 

The  fourth-order  tensor  defined  in  (25)  is  equivalent  to  the  T  matrix  of  quantum  mechanics,  and 
the  nth-order  term  in  its  series  (Bom)  expansion  expresses  the  scattering  process  of  order  n  -  1 . 
Using  (24)  in  (22)  yields  the  exact  solution 

C  =  C°  +  (I  -  (GT))-1  <T) .  (26) 

Because  the  exact  evaluation  of  the  second  term  in  (26)  is  infeasible  (e.g.,  G  must  be  calculated 
for  C(x) ),  we  resort  to  a  self-consistent  approximation:  C°,  which  has  thus  far  been  arbitrary,  is 
chosen  to  equal  C ,  so  that  (26)  is  recovered  if  <T)  =  0 .  Substituting  the  Bom  expansion  in  (25) 
generates  a  sum  in  which  the  nth  term  is  a  contraction  of  the  (n-l)th  outer  product  of  G  against 
the  n-point  moment  of  6C . 

Truncating  the  Bom  series  at  first  order  yields  the  Voigt  average,  C  «  C  .  To^  obtain  a  better 
approximation,  we  have  carried  the  expansion  to  second  order.  Since  C  -  C  ~  0(6C),  the 
second-order  Bom  approximation  takes  the  form  of  an  implicit  integral  equation, 

Q/W  ~  Q/W  +  ^Gpqrs(x)‘Vijpqrski. *) (x ) ,  (27) 

where  Gpqrs  is  the  strain  Green's  function  for  a  transversely  isotropic  medium  with  a  constant 
elasticity  tensor  Cijki,  and  V  ljklpqrs  are  the  components  of  the  eighth-order,  spatially  stationary 
covariance  tensor  V  (x)  =  Ccc(x), 

^ ijklpqr /x)  =  ([Qjd  x’ )  -  Qjki\\Cpqrlx'  +  x)-Cpqrs\  > .  (28) 

^(x)  can  be  evaluated  from  the  Backus  harmonic  decomposition  (12)  by  integrating  two-point 
products  of  the  harmonic  polynomials  against  the  joint  p.d.f.  p(  s'  ,s'  +s).  This  allows  us  to  fully 
exploit  the  transversely  isotropic  stochastic  symmetry  of  the  C  field.  Moreover,  we  can  obtain 
considerable  simplifications  by  rewriting  the  second  term  in  (27)  as  an  integral  over  the 
wavevector  k;  e.g.,  the  3D  Fourier  transform  of  G(x)  depends  only  on  k ,  and  the  spatial 
derivatives  in  (12)  reduce  to  multiplications  by  .  We  are  currently  developing  an  iterative 
method  for  calculating  (27)  using  this  approach. 
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